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Abstract In the frame of a well established lattice gas 
model for granular compaction, we investigate the high 
intensity tapping regime where a pile expands signifi- 
cantly during external excitation. We find that this model 
shows the same general trends as more sophisticated 
models based on molecular dynamic type simulations. In 
particular, a minimum in packing fraction as a function 
of tapping strength is observed in the reversible branch 
of an annealed tapping protocol. 

Keywords Frustrated Lattice Gas • Granular Com- 
paction • Tapping • Reversible Branch 



1 Introduction 

The phenomenology of compaction of granular matter 
under vertical tapping captures the attention of a num- 
ber of investigations and is not exempt from controver- 
sies [l[ . Most studies nowadays focus on the dynamics of 
compaction — the evolution of structural properties as a 
function of the number of taps applied to the sample. In 
spite of being of chief importance, much less work is done 
on the steady state regime achieved by the sample after 
a very large number of taps. At low tapping intensities, 
the relaxation dynamics of these systems is extremely 
slow, which makes the steady state very hard to reach. 
In a pioneering work, Nowak et al. showed that the 
steady state can be achieved by means of a suitable an- 
nealing. Recently, Ribiere et al. 3 argued that the steady 

Paula A. Gago: Departamento de Ffsica, Facultad de Cien- 
cias Exactas, Universidad Nacional de La Plata, Casilla de 
Correo 67, 1900 La Plata, Argentina. 

Nicolas E. Bueno: Laboratorio de Ensenanza de la Ffsica, 
Departamento de Ffsica, Facultad de Ciencias Exactas, Uni- 
versidad Nacional de La Plata, Casilla de Correo 67, 1900 La 
Plata, Argentina. 

L. A. Pugnaloni: Instituto de Ffsica de Lfquidos y Sistemas 
Biologicos (CONICET La Plata - UNLP), Casilla de Correo 
565, 1900 La Plata, Argentina. Tel.: +54-221-4233283, Fax: 
+54-221-4257317, E-mail: luis@iflysib.unlp.edu.ar 



state is indeed obtainable, reproducible and may consti- 
tute a true thermodynamic state for granular systems. 
These experiments show that the packing fraction 4> in 
the steady state is a monotonic decreasing function of the 
tapping intensity: the so-called "reversible branch". A 
number of studies, from experimental @; 3 0] and mod- 
eling [1; S 0; S H; EH EH approaches have confirmed this 
general trend. One exception has been found in a model 
of pentagon packings [Hl;[l3[ where the reversible branch 
presents a monotonic increase of </> as tapping intensity 
is increased. 

Recently, by using a few models of granular depo- 
sition, it has been shown [l4| that <p, in the reversible 
branch, presents a minimum at relatively high tapping 
intensities. This has now been observed in the labora- 
tory [l5[ . Previous experimental studies where unable to 
observe this feature presumably due to the high inten- 
sity taps required. According to the simulation models, 
it may be necessary to transfer enough energy to the 
packing during a tap to induce the system to expand up 
to five times its volume before this minimum can be ob- 
served. This is a rather difficult experiment to carry out 
with a setup designed to explore slow relaxation by gentle 
tapping. In twodimensional systems, the minimum can 
be observed at much lower tapping strength; an effec- 
tive expansion during the tap of less than twice the bed 
height. It is worthmentioning that some of the models 
used in Ref. [lij where the same models used by others 
§ [1] , but focusing in a different region of the parameters 
that control tap intensity. 

The use of annealed tapping to investigate the hys- 
teresis of granular assemblies subjected to tapping was 
also used in the frame of a very simple frustrated lattice 
gas model proposed by Nicodemi et al. 0; [l(| to inves- 
tigate granular compaction. In this work we assess the 
ability of this model to display a reversible branch simi- 
lar to the one observed in experiments and whether the 
density minimum predicted by more sophisticated mod- 
els is also present. With this aim, we explore a range of 
parameter values of the model not yet reported in the 
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literature in order to simulate strong intensity tapping 
conditions. 

The rest of the paper is organized as follows. In Sect. [2] 
we review the main features of the frustrated lattice gas 
model In Sect. [3] we describe the tapping protocol. 
In Sect. S] we present the results. In Sect. [5] we draw the 
conclusions. 



2 Frustrated lattice gas model 

We consider a model first proposed by Nicodemi et al. [9[. 
The model consists of a system of particles that move on 
a square lattice whose bonds are characterized by fixed 
random numbers e^j = ±1. Any site i can be cither 
empty or filled with a particle. Particles are character- 
ized by an internal degree of freedom (spin) Si = ±1 
and are subjected to the constraint that whenever two 
neighboring sites i and j hold particles with spin Si and 
Sj , the following conditions must be satisfied 

CijSiSj = 1. (1) 

This condition introduces an effective frustration. In 
such system there will always be empty sites due to this 
constraint. 

In order to introduce external vibrations and gravity, 
the square lattice is tilted by 45 degree. In this way parti- 
cles diffusing on the lattice will always move up or down 
since there are no neighboring sites at the same vertical 
position. The model is used to carry out a Monte Carlo 
(MC) simulation. In each MC step, we choose a particle 
and a move type (spin flip or particle hopping) at ran- 
dom. A spin flip is accepted with probability one if there 
is no violation of Eq. [T] If a particle attempts a hop, one 
of the four neighboring sites is chosen with equal proba- 
bility as the landing site. If the particle attempts a move 
upward (either up-right or up- left), the landing site is 
empty and the internal degrees of freedom satisfy Eq. [TJ 
we accept the move with probability Pi. If the particle 
attempts to move downward, we accept the move with 
probability P\ (with Pi + P 2 = 1), provided that Eq. [T] 
holds and that the landing site is available. 

In the absence of vibrations, the effect of gravity is 
simulated by imposing Pi = 0. When vibrations are 
switched on, Pi becomes finite. The control parameter is 
the ratio x = PijP\ which describes the intensity of the 
vibration. 

The MC simulations of the model described above are 
done on a tilted lattice with periodic boundary condi- 
tions along the horizontal axis and rigid walls at bottom 
and top. The top wall is always high enough to avoid 
particles to reach it during vibration. After fixing the 
random bond matrix tij, an initial configuration is pre- 
pared by randomly inserting particles of given spin one 
at a time into the box from its top. Each particle falls 
down following the described dynamics with x = 0. After 
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time (2N MC trials) 

Fig. 1 Packing fraction <j> as a function of MC time t during 
the deposition phase in a tap for xq = 1.0 and r = 10 5 . At t = 
10° the vibration is switched off and deposition begins. After 
stabilization, the vibration is switched back on at t m 1.017 x 
10 J . The inset shows <j> along two full taps. The deposition 
phases cover a very narrow slot of time as compared with the 
long duration of the vibration phases in this example. 

a particle has settled, a new one is released. We intro- 
duce a fixed number N of particles. To obtain an initial 
low density configuration, particle spins are not allowed 
to flip in this preparation process. The state prepared 
in this way has a packing fraction of <f> = 0.563 ± 0.007. 
Packing fraction cf> is always measured as the fraction 
of occupied sites in a rectangular region that span the 
system in the horizontal direction and whose height is 
equivalent to the height that the system would have had 
at <f> = 1 (i.e., if no voids where present). The region 
of measurement is vertically centered with the center of 
mass of the system. 



3 Annealed tapping 

In the studies carried out by Nicodemi et al. [l0[ , tapping 
is simulated through different protocols depending on 
the shape of the function x(t) used to control vibrations. 
Here, we use a stepwise function that switches on the 
vibration to a constant value of x over a period of time 
r and then switches it back off until full deposition and 
equilibration is reached. Hence, each tap is simulated 
using 

x(t) = x Q [l - 6(t - t)}, (2) 

with 9(t) the Heaviside step function. 

Full deposition and equilibration is defined as the 
state in which all spins remain in the same site for at 
least 200 unit time even though they may flip. A unit 
time is defined as 2N MC trials. In average, each parti- 
cle has one chance to hop and one chance to flip in this 
time. 
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number of taps + 1 

Fig. 2 Packing fraction if as a function of the number of 
taps applied for xq — 1 and r = 10 (upper curve), 10 4 (lower 
curve), and 1.5 x 10 s (middle curve ) . 

In Fig. [TJ we show the evolution of the packing frac- 
tion 4> during a 'strong' tap with xq = 1.0 and r = 10 5 . 
Whenever the vibration is on [x(t) — xq], 4> (measured in 
a central box, as described in Sect. [5]) fluctuates rapidly 
but continuously decreases in average. Once vibrations 
are turned off [x(t) = 0], <j> rapidly grows as particles 
fall until full deposition is reached and a constant value 
of <j) is obtained. We have run some trial simulations in 
which we wait 1000 unit time (instead of 200) after all 
particle migrations have ceased before deciding that the 
system is fully stable and have found no differences in 
the results. 

Since we are interested in simulating strong taps, 
we have investigated different values of the parameters 
Xq and r in order to achieve an effective expansion of 
the system during vibration. We have found that set- 
ting xq = 1 and varying r allows us to explore a wide 
range of effective expansions. It is worth mentioning that 
Nicodemi et al. [H;[l(| have investigated conditions where 
the slow relaxation of the system were under scrutiny and 
for this reason they were mainly focused on values of Xq 
below 1. We find that for xq = 1 and r < 100 the system 
barely expands during the vibration phase. However, for 
r > 100 an effective expansion of the system is achieved. 
We measure the dimensionless expansion A as the ratio 
between the packing fraction of the system at the time 
we switch off the vibration (t = r) and the packing frac- 
tion after the bed is fully deposited (see Sect. |4]). 

In Fig. [2j we show, as an example, the evolution of 
the packing fraction cf) during tapping with xq = 1 for 
three values of r. These correspond to independent ex- 
periments with an initial low density configuration pre- 
pared as described above. It is clear from the figure that 
low values of t lead to slow relaxation whereas high val- 
ues of t yield a rapid equilibration. 

The annealing protocol is carried out as follows. After 
the initial preparation of the pack, 10 3 taps with xq = 1 
are applied at a low tapping intensity, i.e., with a low 



value of t. Then, the tapping strength is increased by 
increasing the value of r and a new series of 10 3 taps are 
applied to the final configuration reached in the previous 
series. The process is repeated up to very large values of 
r. Then, a decreasing ramp of r values is followed down 
to the initial low tapping intensity. 



4 Results for the high intensity tapping regime 

We use a 30 x 4000 squared lattice where we introduce 
N = 1000 particles. Initially, particles are "dropped" one 
at a time from above (with random horizontal position 
and spin) and they move according to the prescribed dy- 
namics with x = but spin flips are not allowed. When- 
ever a particle is unable to migrate to a neighboring site 
after 16 MC trials the particle is frozen in its site and a 
new particle is released until all particles are deposited. 

The tapping intensity — controlled by r — is increased 
and decreased in discrete steps form r = 1 up to r = 
1.5 x 10 5 and back down to r = 1. At each value of t, a 
total of 1000 taps are applied to the system. The value of 
<p considered for each r corresponds to the average over 
the last 500 taps. The entire annealing is repeated from 
the initial low density configuration for five different ran- 
dom matrices e% j. All values of <j) reported correspond to 
the average over these five independent experiments and 
error bars display the corresponding standard deviation. 
We notice that, although the mean value of 4> has a very 
low dispersion, the "instantaneous" 4> in any given run 
may present larger fluctuations around the mean (see 
Fig. 03 . 

In Fig. [3] we show the packing fraction ^ as a func- 
tion of t for the entire annealing protocol. The increasing 
(irreversible) and decreasing (reversible) branches are in- 
dicated with different symbols. We have also increased 
r from the last high density configuration in a second 
annealing cycle (data not shown). The same values of 
<p obtained during the decreasing ramp are reproduced. 
This confirms that the upper curve shown in Fig. [3] is 
indeed reversible. 

It is clear from Fig. [3] that a minimum in density 
is present at high values of r. Beyond this minimum, (j> 
increases with r rather slowly (notice the log scale in the 
horizontal axis). Due to limitations in CPU time we are 
unable to explore, at present, this regime beyond this 
values of r. However, we can estimate a limiting value of 
(j) for very strong tapping (see below). 

The existence of a minimum in the density-tapping 
intensity curve has been reported very recently |14| . In 
Ref. [14(, the effect is explained in terms of the forma- 
tion of arches. At low tapping intensities, the introduc- 
tion of free volume is very limited (small expansion) . An 
increase in effective expansion promotes the formation 
of larger arches (and hence larger voids) since particles 
need some space in order to arrange in these coopera- 
tive structures. However, if expansion during tapping is 
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Fig. 3 Packing fraction (f> as a function of r for xo = 1 
along an annealing protocol. The circles correspond to the 
"non-reversible branch" of initial increasing tapping intensi- 
ties (increasing r). The squares correspond to the "reversible 
branch" of decreasing tapping intensities (decreasing r). The 
system is tapped 10 3 times at each value of r. The last 
500 configurations are averaged to estimate the value of <j> 
reached. The entire simulation has been repeated five times 
for different bond random matrices Eij. The reported data 
corresponds to the average — and the error bars to the stan- 
dard deviation — over these five repetitions 
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Fig. 4 Packing fraction <f> as a function of the effective di- 
mensionless expansion A induced by tapping. For each value 
of r the corresponding expansion is calculated as the ratio be- 
tween the density at the time when the vibration is switched 
off and the density after full deposition. Symbols as in Fig. [3] 
We show results taken from Ref. 14] for a molecular dynamics 
of dissipative soft disks with static friction (down triangles), a 
pseudo dynamic model of inelastic disks (right triangles) and 
an off lattice MC of hard disks (up triangles). The horizontal 
dash-dotted line corresponds to the limit of very strong taps 
(r —5- oo ) in the lattice model (see text). 



too large, falling particles have less chances to meet each 
other at the time of reaching their stable positions in 
order to cooperate and form arches. 

In the present model arches are not clearly defined. 
A definition of which particles support a given parti- 
cle is required to identify arches [TJ [Hj]. However, the 
frustration introduced by the bond matrix eij induces 
the particles to interlock in loops that leave voids in the 
structure. These loops are the analogous to arches in 
more realistic models. 

In order to compare our results with previous mod- 
els, we plot in Fig. 0] the packing fraction as a function 
of the effective dimensionless expansion A induced by 
the corresponding value of r. Results from Ref. [l4| for 
twodimensional models are also reproduced. We observe 
that all models present a minimum even though the po- 
sition varies in the range 1.10 < A < 2.10. The posi- 
tion of the minimum found in 3D models seems to be in 
the range 3.0 < A < 5.0 The frustrated lattice gas 
model yields rather low packing fractions when compared 
with more realistic models. However, it is clear that the 
main qualitative features of the reversible branch of the 
annealed tapping curve are very well captured by the 
lattice model. 

Since the increase in <fi at large values of r seems to be 
due to the particles depositing without much chances to 
interact as they reach the surface of the growing stable 
pile, we can simulate this in a extreme condition. Instead 
of expanding the system using very long r values, we now 
fill the box dropping particles one at time so that they 
reach the pile individually. This is a process similar to 
the one used in generating the very first initial condition. 



However, we now allow the particles to flip their spin as 
they fall. Moreover, the already deposited particles are 
not freezed but allowed to diffuse and flip while the pile 
is being built up. We expect this filling process to mimic 
the limiting case of an expansion obtained with r — > oo. 
The packing fraction so obtained is shown in Fig. @] We 
argue that this should be the limiting value to which <f> 
would approach if we further increased r. 



5 Conclusions 

We have shown that a well known frustrated lattice gas 
model of granular compaction — initially designed to study 
slow dynamics under gentle tapping — is able to repro- 
duce the qualitative features of more complex models in 
which strong tapping conditions are simulated. In par- 
ticular, we observe the existence of a minimum in the re- 
versible branch of the annealed tapping. It is somewhat 
surprising that such simple model can display nonmono- 
tonic behavior. This makes the model a suitable working 
platform to study this newly discovered feature of gran- 
ular compaction. It would be of special interest to inves- 
tigate if the steady states of equal mean packing fraction 
produced at different tapping intensities (at both sides 
of the minimum) can be distinguished through the size 
of the volume fluctuations, an order parameter, or the 
distribution of characteristic structures such as voids or 
loops. This model presents the advantage that is of sim- 
ple implementation and that simulation are much less 
CPU time demanding than more realistic models which 
show the same general features. 
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In Ref. |14j . the existence of a minimum in the re- 
versible branch of granular compaction is explained in 
terms of arching. Since this lattice model does not con- 
sider proper contact dynamics nor realistic volume ex- 
clusion, it is clear that the presence of the minimum is 
due to an underlying phenomenon that controls both, 
arching and packing fraction. 

At low tapping intensities, the high packing fractions 
achieved are due to repeated tapping that allows the 
system to relax in search for progressively more com- 
pact (lower potential energy) structures. This is achieved 
through the local rearrangements promoted from tap to 
tap. However, at very high tapping intensities, the sys- 
tem achieves the steady state in a single tap (see Fig. [2]) 
since the full pack is arranged de novo after a strong tap 
in a global relaxation event. In this last case, in order 
to achieve higher packing fractions, the particles need to 
follow a deposition process that allows them enough time 
to search for the lowest potential energy configuration in 
a single tap. This is best promoted if the bed is expanded 
significantly during the tap so that the longer deposition 
times required and the larger free volume introduced give 
greater chances for the particles to find lower positions. 
In summary, the existence of the minimum in the re- 
versible branch of granular compaction seems to be ulti- 
mately related to a competition between the global and 
local relaxation promoted at different scales of tapping 
strength @. 
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